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Abstract.  Sharp  fronts  with  rapid  changes  in  fluid  saturations  over  short 
distance  and  time  scales  often  exist  in  multiphase  flow  in  subsurface  systems.  Such 
highly  nonlinear  problems  are  notoriously  difficult  to  solve,  and  standard  solution 
approaches  are  often  inefficient  and  unreliable.  We  summarize  four  existing  and  one 
new  transformation  method  (IT2)  for  solving  Richards’  equation  within  a  common 
framework  and  compare  performance  for  a  wide  range  of  media  properties.  We  show  that: 
transformation  methods  can  significantly  improve  solution  efficiency  and  robustness 
compared  to  standard  solution  approaches;  transformation  parameters  depend  upon 
auxiliary  conditions,  media  properties,  and  spatial  and  temporal  discretization;  and  IT2 
compares  favorably  with  existing  transforms. 
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1.  Introduction 

Modeling  flow  of  multiple  fluids  in  subsurface  systems  is  an  important  water 
resources  problem  for  which  many  unresolved  questions  remain  [ Miller  et  al.,  1997]. 
Resolving  sharp  fronts  of  a  dependent  variable  as  a  function  of  space  and  time  is  a 
particularly  challenging  aspect  of  this  class  of  problems.  These  fronts  can  develop 
when  a  more  viscous  wetting  phase  displaces  a  less  viscous  non- wetting  phase,  as  when 
water  infiltrates  a  porous  media  that  is  initially  relatively  dry,  resulting  in  sharp  fronts 
in  volumetric  fractions  of  the  fluid  phases  for  certain  media  properties  and  auxiliary 
conditions. 

While  such  sharp  fronts  result  from  many  possible  combinations  of  conservation 
laws  and  constitutive  relations  of  interest  in  multiphase  flow  problems,  Richards’ 
equation  (RE)  with  common  constitutive  relations  is  a  simple  case  of  substantial 
practical  importance  where  this  difficult  class  of  problems  can  arise  [Miller  and  Kelley, 
1994;  Tocci  et  al.,  1997].  Because  of  this,  RE  is  a  good  test  problem,  although  it  is 
understood  that  significant  theoretical  questions  remain  unresolved  about  its  adequacy 
for  describing  unsaturated  flow  [ Gray  and  Hassanizadeh,  1991a;  Gray  and  Hassanizadeh, 
1991b;  Miller  et  al.,  1997]. 

Sharp  volumetric  fluid-phase  fraction  fronts  result  in  part  from  nonlinearities  in 
constitutive  relations  that  describe  the  interdependence  of  fluid  pressures,  saturations, 
and  conductivities.  These  nonlinearities  are  present  in  discretized  forms  of  the  governing 
equations  that  are  solved.  Accurate  resolution  of  these  sharp  fronts  can  require  fine 
discretizations  in  space  and  time,  and  can  lead  to  systems  of  nonlinear  algebraic 
equations  that  are  difficult  to  solve.  Choosing  discretization  or  nonlinear  solution 
methods  inappropriately  can  lead  to  inefficient  or  unreliable  approximate  solution 
approaches. 

Several  approaches  to  these  problems  can  be  grouped  together  under  the  rubric  of 
transformation  methods.  Transformation  methods  seek  to  reduce  the  sharpness  of  a 
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front  in  a  problem  through  the  identification  and  application  of  an  appropriate  change 
of  variable  applied  to  the  dependent  variable.  The  original  problem’s  solution  may  then 
be  retrieved  by  applying  an  inverse  transformation. 

While  several  transformation  methods  have  been  applied  to  RE  with  some  success, 
we  lack  detailed  investigations,  comparisons,  and  guidance  for  the  applying  these 
approaches  generally.  Further,  some  of  the  transformation  methods  applied  to  date 
include  a  parameter(s)  that  must  be  specified,  and  the  effect  of  the  parameter’s  value 
on  solution  efficiency  and  robustness  is  often  unknown.  These  difficulties  provide  a 
barrier  to  the  routine  application  of  transformation  methods  for  solving  RE  and  other 
multiphase  flow  problems. 

The  goal  of  this  work  is  to  advance  the  current  understanding  of  transformation 
approaches  for  solving  Richards’  equation.  The  specific  objectives  of  this  work  are  (1) 
to  develop  a  framework  to  summarize  transformation  methods  for  solving  RE;  (2)  to 
identify  key  components  of  a  solution  approach  for  RE  that  may  affect  the  efficiency 
and  robustness  of  the  resulting  approximation;  (3)  to  compare  existing  transformation 
methods  for  a  range  of  media  and  auxiliary  conditions;  (4)  to  investigate  the  importance 
of  parameter  selection  on  the  efficiency  of  parameterized  transformation  methods;  (5)  to 
identify  and  evaluate  a  new  transform  for  a  wide  range  of  test  problems;  and  (6)  to  give 
guidance  for  the  application  of  transformation  methods  to  solve  RE. 

2.  Background 

Transformation  methods  for  solving  RE  have  existed  for  several  decades  [Rubin, 
1968;  Raats  and  Gardner,  1974;  Haverkamp  et  al.,  1977;  Baca  et  al.,  1978;  Vauclin  et 
al,  1979].  The  general  objective  of  transformation  methods  is  to  overcome  inefficiencies 
in  the  numerical  solution  process — inefficiencies  caused  by  the  strong  nonlinearity  of  the 
media  hydraulic  properties  as  functions  of  pressure,  particularly  in  the  case  of  infiltration 
into  a  media  that  is  initially  relatively  dry.  These  types  of  infiltration  problems  give  rise 
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to  very  high  water  pressure  gradients  near  the  wetting  front  and  lead  to  computationally 
inefficient  numerical  solutions  when  using  standard  techniques.  Prohibitively  small 
time  step  sizes  or  a  large  number  of  nonlinear  iterations  are  often  required  for  such 
problems.  Combining  transform  methods  with  iterative,  implicit,  mass-conserving 
numerical  methods  results  in  more  efficient  solutions.  Current  transformations  include 
integral  [ Haverkamp  et  al. ,  1977],  hyperbolic  function  [i?oss,  1990],  variable  switching 
[. Kirkland  et  al.,  1992;  Forsyth  et  al,  1995],  and  rational  function  transformations  [ Pan 
and  Wierenga,  1995]. 

Early  attempts  at  transformation  methods  used  an  integral  transform  commonly 
referred  to  as  the  Kirchoff  integral  transformation  (IT1)  [Rubin,  1968;  Raats 
and  Gardner,  1974;  Vauclin  et  al.,  1979;  Haverkamp  et  al.,  1977;  Redinger  et  al., 
1984;  Campbell,  1985].  This  transformation  directly  reduces  the  nonlinearity  of  the 
conductivity  terms  in  RE  and,  as  a  result,  is  effective  in  reducing  the  number  of 
nonlinear  iterations  required  for  a  solution,  under  certain  conditions.  However,  IT1 
depends  on  media  hydraulic  properties  and  will  therefore  vary  spatially  if  the  hydraulic 
properties  do  so.  Thus,  simple  application  of  IT1  is  restricted  to  homogeneous  media. 
Corrections  can  be  made  to  adapt  IT1  to  layered  and  gradational  media  by  adding  a  flux 
balancing  correction  [ Ross  and  Bristow,  1990].  IT1  is  also  more  complex  to  implement, 
since  an  analytic  function  of  the  inverse  is  generally  not  available. 

Solving  the  water-content-based  form  of  RE  can  result  in  significantly  improved 
performance  compared  to  the  traditional  pressure-based  methods  [Hills  et  al.,  1989]. 
This  is  due  to  the  fact  that  the  media  hydraulic  functions  are  less  nonlinear  when 
expressed  in  terms  of  water  content  than  when  expressed  in  terms  of  pressure, 
particularly  when  modeling  infiltration  into  a  relatively  dry  media.  Thus  numerical 
solutions  using  the  water-content-based  RE  generally  require  fewer  nonlinear  iterations. 
The  limitation  of  this  approach  is  that  the  water-content-based  RE  cannot  be  used  to 
solve  infiltration  problems  involving  saturated  regions.  A  0-based  transform  (THT)  is 
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an  attempt  to  retain  the  advantages  of  water-content-based  methods  while  remaining 
applicable  to  media  containing  saturated  or  near-saturated  regions.  THT’s  include 
the  affine  transformation  of  9  approach  [Kirkland  et  al,  1992]  as  well  as  the  primary 
variable  switching  technique  [ Forsyth  et  al. ,  1995].  The  THT  has  the  characteristics 
of  water  content  when  the  media  is  unsaturated  and  of  pressure  when  the  media  is  at 
or  near  saturation.  This  approach  has  shown  roughly  the  same  reduction  in  nonlinear 
iterations  as  in  the  water-content-based  solution.  Since  the  THT  is  defined  in  terms  of 
volumetric  water  fraction,  it  will  vary  spatially  if  the  media  type  does  so.  As  is  the  case 
with  IT1,  simple  application  of  THT  is  restricted  to  homogeneous  media. 

Because  direct  application  of  IT1  and  THT  is  restricted  to  homogeneous  media, 
an  alternate  class  of  transforms  was  developed  that  would  be  directly  applicable  to 
heterogeneous  media.  These  transforms  are  defined  strictly  in  terms  of  ip  and  arbitrary 
parameter(s).  Since  ip  is  continuous  across  the  boundaries  between  different  media 
types,  these  transform  functions  will  be  continuous  variables  in  heterogeneous  media, 
provided  that  the  parameters  remain  constant  across  media  types.  The  first  such 
transform  dates  back  to  a  simple  log  transform  [ Baca  et  al.,  1978].  This  was  effective, 
but  a  more  general  and  efficient  transform,  defined  in  terms  of  the  hyperbolic  sine 
function  (HST),  was  subsequently  introduced  [Ross,  1990].  HST  reduces  computational 
expense,  but  it  introduces  two  arbitrary  parameters.  It  is  recommended  [Ross,  1990] 
that  one  of  the  parameters  be  fixed  at  a  constant  percentage  of  the  other  to  reduce  HST 
to  one  arbitrary  parameter  [Ross,  1990].  However,  to  optimize  the  HST  completely  for 
a  given  problem,  the  two  parameters  must  remain  arbitrary.  Since  the  introduction  of 
HST,  another  of  this  class  of  transforms  has  been  proposed.  This  transform  is  defined  in 
terms  of  a  rational  function  (RFT)  of  ip  [Pan  and  Wierenga,  1995].  The  RFT  provides 
performance  improvements  similar  to  HST,  but  introduces  only  one  arbitrary  parameter. 

Because  all  of  the  transforms  except  IT1  involve  arbitrary  parameters,  selecting 
parameter  values  is  important  to  determine  the  efficiency  of  a  particular  transformation. 
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The  literature  to  date  lacks  rigorous  optimization  of  transform  parameters  over  a 
wide  range  of  test  problems.  The  literature  suggests  that  selection  of  HST  [Ross, 
1990]  and  RFT  parameters  [Pan  and  Wierenga,  1995]  for  optimal  or  near-optimal 
performance  is  independent  of  media  properties,  but  this  claim  must  be  verified  through 
formal  optimization  and  sensitivity  analysis.  Such  analysis  is  critical  in  determining 
the  efficiency  and  robustness  of  a  given  transformation  but  has  not  been  reported  in 
the  literature  to  date.  If  the  data  suggest  a  relationship  between  effective  transform 
parameters  and  media  properties  for  the  HST  and  RFT  type  transformations,  then  the 
advantage  that  these  transforms  are  applicable  in  their  current  form  to  heterogeneous 
media  conditions  becomes  less  certain. 


3.  Formulation 


3.1.  ID  Richards’  Equation 


We  consider  one-dimensional  (ID)  infiltration  in  this  work,  beginning  with  a 
general,  p  version  of  the  ID  RE.  For  the  case  in  which  fluid  compressibility  is  included 
for  a  vertical  system,  the  p  version  of  the  common  mixed  form  equation  is  given  by 


dtp  dp 

SsSa  ^  dp  dt  + 


d6  (p) 
dt 


d_ 

dz 


K(p ) 


(  dip  dp  \ 

J 


(1) 


where  Ss  is  the  specific  storage  coefficient,  which  accounts  for  fluid  compressibility;  Sa 
is  saturation  of  the  aqueous  phase;  p  —  p{ip)  is  a  general  transformation  function;  ip  is 
the  pressure  head;  t  is  time;  9  is  the  volumetric  water  fraction  of  the  aqueous  phase;  z  is 
the  vertical  spatial  dimension;  and  K  is  the  hydraulic  conductivity. 

We  consider  problems  with  auxiliary  conditions  of  the  form 


p(z,t  =  o) 

=  Po  (z)  =  p  (ipo  (z)) 

(2) 

p  (z  =  0,  t  >  0) 

=  Pi  =p  (ipl) 

(3) 

p  (z  —  Z,  t  >  0) 

=  P2=P  M 

(4) 
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where  Z  is  the  length  of  the  domain,  ip0  may  be  a  function  of  space,  and  ipi  and  xp2  are 
constants.  The  auxiliary  conditions  p0l  pi,  and  p2  are  found  by  the  appropriate  change 
of  variables  of  the  initial  and  boundary  conditions  given  in  terms  of  ip. 


3.2.  Constitutive  Relations 


Solution  of  RE  requires  constitutive  relations  to  describe  the  interdependence 
among  pressure,  saturation,  and  hydraulic  conductivity.  For  the  purpose  of  this  work, 
we  use  the  well-known  van  Genuchten  (VG)  pressure-saturation  relationship  [ van 
Genuchten,  1980],  given  by 

9a  W  -9r  _  j  (1  +  <  0 

1,  ip  >  0 


Se  W  = 


9 ,  -  9r 


(5) 


where  mv  —  1  —  l/nv]  Se  is  the  effective  saturation;  9r  is  the  residual  volumetric  water 
content;  9S  is  the  saturated  volumetric  water  content;  and  av  and  nv  are  experimentally 
determined  coefficients  in  the  VG  p-S  model,  which  are  related  to  the  mean  pore  size 
and  the  uniformity  of  the  pore  size  distribution,  respectively. 

The  saturation-permeability  relation  for  the  aqueous,  or  wetting  phase,  is  described 
using  Mualem’s  model  [ Mualem ,  1976] 

K  (Se)  =  KsSlJ2  [l  -  (l  -  (6) 

where  Ks  is  the  water-saturated  hydraulic  conductivity,  and  Se  =  Se(xp)  from  (5). 


3.3.  Transform  Definitions 

Starting  from  the  general  p  version  of  the  ID  RE  (1),  the  more  commonly  used 
mixed,  ^-based,  and  0-based  (water-content-based)  forms  of  RE  can  be  retrieved.  The 
basic  idea  behind  the  transformation  approach  is  to  define  a  function  p(ip)  that  will 
result  in  a  more  efficient  and  robust  solution  to  the  governing  equation,  (1).  Several 
transformations  have  been  developed  to  date,  and  the  definitions  of  p(ip)  for  each  are 
listed  in  Table  1.  Arbitrary  parameters  are  represented  by  (3  and  (i\. 


Table  1 
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THT  and  IT1  linearize  either  the  temporal  or  spatial  derivative  term  in  RE, 
respectively.  While  it  is  not  possible  to  linearize  both  the  temporal  and  spatial  terms  in 
RE  completely,  the  efficiency  and  accuracy  of  the  numerical  solution  can  be  significantly 
enhanced  by  choosing  a  transform  that  achieves  the  proper  balance  between  linearization 
of  temporal  and  spatial  terms.  RFT  and  HST  can  partially  linearize  both  temporal  and 
spatial  terms;  they  are  defined  explicitly  in  terms  of  xp  instead  of  K  or  9.  This  allows 
for  a  straightforward  application  for  heterogeneous  media,  but  such  applications  will 
be  effective  only  if  the  arbitrary  parameter  is  not  dependent  upon  and  sensitive  to  the 
given  media  parameters. 

As  an  alternative  to  extant  transforms,  we  tested  several  alternative  transforms 
in  an  attempt  to  find  one  that  would  effectively  balance  the  linearization  of  temporal 
and  spatial  terms  and  would  not  be  as  sensitive  to  changes  in  media  parameters.  We 
evaluated  higher-order  rational  functions,  exponential-based  functions,  and  integral 
functions.  We  looked  at  transforms  defined  explicitly  in  terms  of  ip  as  well  as  those 
that  included  media  hydraulic  functions  K  and  9  in  their  definitions.  Based  upon  our 
evaluation,  we  have  found  an  integral  transform  that  is  effective  in  terms  of  efficiency 
and  robustness.  This  transform,  IT2,  is  defined  in  terms  of  9  and  the  integral  of  K. 

Note  that  all  of  the  transform  functions  are  C 1  continuous  functions  of  xp  and  have 
non-zero  derivatives  for  all  values  of  xp]  i.e. ,  dp/dxp  ^  0.  These  conditions  are  necessary 
for  accurate  and  efficient  solution  of  (1).  Note  also  that  for  the  case  of  (3  =  0,  IT2 
reduces  to  IT1.  Thus,  IT1  is  a  subset  of  IT2. 

4.  APPROACHES 

4.1.  Spatial  and  Temporal  Approximation 

The  numerical  solution  of  (1)  begins  with  a  discrete  formulation,  for  which  we  use  a 
spatially-centered,  fully-implicit  finite  difference  approximation.  For  the  spatial  domain 
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[0,Z],  a  uniform  grid  is  defined  with  nn  —  1  intervals  [zj,  ^+1]^1_1,  of  length  A z,  with 
A;?  =  Z/(nn  —  1),  and  Z{  —  [i  —  1)A^  for  1  <  *  <  nn.  (1)  is  then  approximated  for 
1  <  i  <  nn  as 


where  nn  is  the  number  of  spatial  nodes  in  the  solution,  p%  is  the  approximation  to  the 
solution  pfa),  l  indicates  time  level,  and  K^y2  are  interblock  conductivities  estimated 
from  known  values  of  K{p\ |f )«,  K{p\+l),  and  K{p\%\). 

Of  the  several  ways  to  estimate  interblock  conductivities  [ Haverkamp  and  Vauclin, 
1979;  Zaidel  and  Russo,  1992;  Warrick,  1991],  we  investigated  four  approaches.  The 
first  approach  we  used  was  the  arithmetic  mean  technique  (KAM): 

Ki± i/2  =  (Ki  +  Ki± i)  /2  (8) 

KAM  has  been  used  routinely  [ Haverkamp  and  Vauclin,  1979;  Warrick,  1991;  Zaidel 
and  Russo,  1992]  and  is  simple  and  inexpensive  to  compute. 

Due  to  the  nonlinearity  of  ip,  the  interblock  conductivity  estimation  technique 
defined  by  K[(ipi  +  tpi± i)/2]  is  ineffective.  Nonlinear  solvers  will  often  fail  at  small 
simulation  times  when  using  this  approach.  But  since  the  transformed  variable  p  will  be 
somewhat  smoother  in  terms  of  the  spatial  dimension,  the  second  approach  we  used  was 
the  transformed  analog  of  this  approach,  i.e.,  conductivity  evaluation  at  the  arithmetic 
mean  of  the  transformed  variable  p  (KAMP).  KAMP  is  defined  as 

Ki±  1/2  =  K  [( pi  +  pi± i)  /2]  (9) 

Because  K  varies  in  space  as  a  function  of  ip,  an  integral  representation  of  mean 
interblock  values  can  be  computed  [  Warrick ,  1991;  Zaidel  and  Russo,  1992],  The  third 
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approach  we  used  is  an  integral  representation  in  transformed  space  (KINT): 


K%±  1/2  —  < 


rmax{pi,pi±1}  , 


I  .  -  ~Y(K^-)dp 

I'nautPj.PiiUWU  ’ 

Jmin{pi,Pi±i  }  \  dp  J 

l  Ku 


if  Pi  ±  Pi± i 
if  pi  =  pi± i 


(10) 


This  approach  has  not  been  widely  used  for  the  solution  of  RE  because  of  the  apparent 
expense  involved  in  computing  the  integrals.  This  approach  has  not  been  used  in  a 
transform  solution  to  RE  to  our  knowledge. 

The  fourth  approach  we  used  for  estimating  interblock  conductivities  is  termed  the 
arithmetic  mean  saturation  (KAMS)  [ Zaidel  and  Russo,  1992]  and  is  computed  by 


if, ±1/2  =  K  [(s,  +  S,±1)  /2] 


(11) 


This  technique  is  easy  to  implement  but  if  spline  approximations  are  used,  this  approach 
is  most  efficiently  implemented  by  splining  in  terms  of  Se.  which  itself  is  splined  in  terms 
of  p.  This  results  in  some  added  expense.  Spline  issues  are  discussed  in  more  detail 
below. 


4.2.  Nonlinear  Solution  Methods 

To  find  an  approximating  solution  {pp,  1  <  i  <  nn}  to  (7)  requires  solving  a  system 
of  nonlinear  equations.  For  the  purpose  of  this  work,  we  use  either  modified  Picard 
iteration  (MPI)  [ Celia  et  al,  1990]  or  full  Newton  iteration  (NI).  The  MPI  technique 
produces  a  simple  computational  algorithm  that  is  mass  conserving  for  numerical 
approximations  that  preserve  spatial  symmetry. 

The  NI  system  is  formed  from  (7)  by  defining 


(12) 
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and  then  solving  for  f(p)  =  0  by  Newton  iteration 

[J]{Ap}  =  -{/}  (13) 

where  [J]  is  the  Jacobian.  NI  applied  to  a  solution  algorithm  that  uses  fixed  time  steps 
is  not  an  effective  strategy,  since  reliable  convergence  often  requires  very  small  time 
steps.  We  overcome  this  problem  by  using  a  quadratic/cubic  line  search  technique 
[Dennis  and  Schnabel,  1996]  in  concert  with  NI,  which  we  term  a  Nl-line  search  (NILS) 
approach.  NILS  is  much  more  robust  than  NI. 

We  terminated  the  nonlinear  solvers  based  on  an  absolute  error  measure  of  p.  Thus, 
if  the  ranges  of  the  various  transforms  are  not  equal,  those  with  smaller  ranges  will  tend 
to  converge  sooner,  though  often  with  increased  error.  To  avoid  this  potential  problem, 
the  transforms  were  normalized  so  that  the  ranges  of  each  are  equivalent;  i.e.,  each 
of  the  transform  functions  maps  the  domain  [ipmin,ipm\  onto  the  range  [0.pm\,  where 
ipm  =  mini 0,^max],  and  pm  —  1pm  Ipmin- 

4.3.  Efficiency  Considerations  and  Evaluation 

We  define  efficiency  as  the  computational  effort  required  to  achieve  a  specified 
accuracy.  The  evaluation  of  constitutive  relations  requires  a  significant  amount  of 
computational  effort  using  the  standard  approach  for  solving  RE.  Therefore  these 
relations  are  often  evaluated,  the  values  stored  in  tables,  and  intermediate  values 
determined  by  linear,  cubic,  or  Hermite  spline  interpolation.  This  procedure  results  in 
a  significant  savings  in  computational  effort  compared  to  direct  function  evaluations 
without  a  significant  change  in  the  accuracy  of  the  solution.  Based  upon  our  previous 
work  and  some  additional  screening  for  transformed  solutions,  we  used  Hermite  splines 
in  this  work,  which  are  described  in  detail  elsewhere  [Miller  et  al.,\. 

Table  interpolation  reduces  computational  effort  even  more  for  transformation 
approaches  than  for  standard  approaches  for  solving  RE.  This  is  especially  so  for  the 
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more  complex  transforms  such  as  IT2,  which  would  require  the  determination  of  roots 
within  a  numerical  integration  procedure.  Table  interpolation  greatly  reduces  the 
computational  effort  per  iteration,  although  this  effort  is  still  greater  than  the  effort 
required  for  a  table  interpolation  approach  to  an  untransformed  solution.  The  difference 
in  effort  between  a  transformed  and  an  untransformed  approach  depends  upon  the 
nonlinear  solution  scheme  used,  although  this  difference  is  usually  less  than  50%  per 
iteration.  Clearly  then,  transformation  methods  must  require  fewer  iterations  than 
untransformed  solutions  in  order  to  be  more  efficient. 

Error  vs.  work  plots  are  used  to  compare  the  various  transforms  in  terms  of 
efficiency,  robustness,  and  parameter  sensitivity.  These  plots  not  only  show  the  overall 
effectiveness  of  a  given  transform  in  terms  of  speed  and  accuracy  over  a  wide  range  of 
media  conditions,  they  also  serve  to  illustrate  optimal  parameter  range  and  parameter 
sensitivity.  To  compare  the  transformation  approaches  to  the  traditional  untransformed 
approach,  the  cost  per  iteration  used  to  calculate  work  as  a  function  of  nonlinear 
iterations  is  adjusted  to  account  for  the  additional  effort  required  in  the  transformed 
solution  approach. 

For  methods  based  upon  the  MPI  approach,  the  work  primarily  concerns  forming 
the  coefficient  matrix  and  right  hand  side  vector,  and  solving  the  linear  systems  of 
equations.  This  observation  allows  for  a  simple,  straightforward  measure  of  work  that 
requires  relative  weights  for  the  two  procedures  and  integer  counts  for  each  of  the 
procedures,  such  as 

Wp  =  wcnc  +  wtni  (14) 

where  Wp  is  a  work  measure  for  MPI  methods,  wc  is  a  weighting  factor  for  formation 
of  the  coefficient  matrix  and  right  hand  side  vector,  which  are  typically  done  at  the 
same  time,  wi  is  a  weighting  factor  for  solution  of  the  linear  system  of  equations,  nc 
is  the  number  of  coefficient  matrix  formation  calls,  and  rq  is  the  number  of  linear 
solutions  performed.  As  reported  in  previous  work  [Tocci  et  al,  1997],  estimates  of 
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the  weighting  coefficients  based  on  detailed  profiling  results  for  the  untransformed  RE 
using  MPI  solver,  KINT  permeability  approximations,  and  Hermite  spline  interpolation 
are  (wc)ut  =  0.530  and  (wi)ut  =  0.181.  We  performed  similar  profiling  tests  using  the 
transformation  approach  and,  based  on  these  results,  our  estimates  for  the  solution  of 
the  transformed  RE  using  MPI,  KINT,  and  Hermite  splines  are  (wc)tr  —  0.743  and 
(wi)tr  —  0.181.  Thus  for  this  case,  the  transformation  approach  requires  about  30% 
more  effort  per  iteration  than  the  untransformed  solution. 

Error  was  evaluated  by  comparison  to  a  dense-grid  solution.  This  error,  referred  to 
as  dense-grid  error,  is  defined  by 
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where  k  is  the  norm  measure  and  i/i  is  an  accurate  approximation  of  the  true  solution 
based  on  a  dense  spatial  grid,  k  =  l^  k  =  2,  and  k  —  oo  were  considered  in  this  work 
and  termed  L2l  and  L <*,  error  norms,  respectively.  The  dense  grid  solutions  were 
generated  using  the  MPI  solver  with  temporal  and  spatial  grid  sizes  equal  to  1/32  of  the 
standard  sizes  listed  in  Tables  2  and  3. 


4.4.  Parameter  Optimization 

The  parameters  for  any  of  the  transformations  can  be  optimized  for  some 
performance-based  objective  function  such  as  amount  of  work  required  or  dense-grid 
error.  In  this  work,  we  use  the  objective  function 
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where  /3  is  the  arbitrary  transform  parameter,  Wp  is  the  required  work  as  defined  by 
Equation  (14),  and  ||  eD  ||i  is  the  dense-grid  error  as  defined  by  Equation  (15). 

For  the  parameter  optimization  of  the  work,  we  used  the  nonlinear  optimization 
package  IFFCO  [Gilmore  and  Kelley,  1993].  IFFCO  is  a  projected  quasi-Newton 


15 


algorithm  that  uses  a  decreasing  sequence  of  finite  difference  steps  (scales)  to 
approximate  the  gradient.  It  uses  an  approximation  to  the  Hessian  and  a  line  search 
algorithm  that  gives  the  code  global  convergence  capabilities. 

4.5.  Test  Problems 

We  compared  the  transformation  approaches  for  solving  RE  to  traditional  solution 
methods  using  eight  sets  of  test  conditions,  which  are  summarized  in  Tables  2  and  3. 
Four  of  these  conditions  (Problems  A-D)  have  been  used  previously  in  the  literature 
as  test  problems  [Miller  and  Kelley,  1994]  [Celia  et  al,  1990;  Rathfelder  and  Abriola, 
1994]  [Forsyth  et  al.,  1995].  Problems  EH  represent  various  soil  textural  groups  sand, 
loamy  sand,  loam,  and  clay  loam,  respectively,  according  to  the  USDA  classification 
[Soil  Conservation  Service,  1975]  as  estimated  by  Carsel  and  Parrish  [Carsel  and 
Parrish,  1988]  from  analyses  of  a  large  number  of  soils. 

With  the  exception  of  Problems  C  and  D,  each  set  of  simulation  conditions  yields  a 
difficult  sharp-front  problem  with  relatively  dry  initial  conditions.  Problems  C  and  D 
were  considered  because  they  offer  an  excellent  benchmark  for  comparing  our  results 
to  recent  research  performed  using  state-of-the-art  methods.  However,  Problems  C 
and  D  are  substantially  easier  than  the  remaining  problems  because  the  domain  is 
much  smaller,  the  initial  conditions  are  less  severe,  and  (for  Problem  C)  fully  saturated 
conditions  do  not  develop. 

To  illustrate  the  effect  of  the  transformation  approach  on  resulting  solution  profiles, 
sample  results  are  plotted  for  Problem  A.  Figure  1  shows  the  solution  in  terms  of  the 
untransformed  pressure  head  ip  at  various  simulation  times.  Figure  2  shows  the  solution 
in  terms  of  the  transformed  variable  p  at  various  simulation  times,  using  the  optimized 
IT2  transform.  Clearly,  the  infiltration  profile  is  smoother  when  plotted  in  terms  of  the 
transformed  variable  p. 


Tables  2  and 


Figure  1 
Figure  2 
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5.  Results  and  Discussion 

The  number  of  transforms  (four),  interblock  conductivity  methods  (four),  nonlinear 
solution  methods  (two),  test  problems  (eight),  and  other  variables  of  concern,  such 
as  spatial  and  temporal  discretization  combine  to  yield  a  large  number  of  possible 
combinations.  While  several  thousand  simulations  were  performed  in  this  work,  a 
complete  analysis  of  each  of  these  variables  was  beyond  the  scope  of  this  effort.  In  the 
sections  that  follow  we  report  on  results  of:  (1)  baseline  comparisons  for  interblock 
conductivity  and  nonlinear  solution  methods;  (2)  performance  comparisons  for  all 
transforms  and  test  problems  using  a  single  nonlinear  solver  and  interblock  conductivity 
method;  and  (3)  the  sensitivity  of  the  transform  parameter,  [J,  to  spatial  and  temporal 
discretization  levels. 

5.1.  Baseline  Comparisons 

For  each  of  the  eight  test  problems,  baseline  comparisons  were  made  for  simulations 
incorporating  all  possible  combinations  of  transformation,  interblock  conductivity 
estimation,  and  nonlinear  solver.  All  of  the  simulations  were  made  using  fixed  time 
steps  and  one  arbitrary  transformation  parameter.  For  the  case  of  the  HST  this  was 
accomplished  by  setting  fli  =  0.1  (3  [i?oss,  1990].  Parameter  optimization  was  performed 
for  all  transforms  on  each  test  problem  using  IFFCO.  These  results  were  used  to  identify 
the  most  promising  combinations  of  interblock  conductivity  approach,  nonlinear  solver, 
and  transform. 

While  these  results  are  not  reported  in  detail,  we  draw  the  following  general 
conclusions  from  our  baseline  runs:  (1)  the  integral  conductivity  estimation  method 
consistently  demonstrated  a  greater  reduction  in  dense-grid  error  than  the  other  three 
interblock  conductivity  approaches  with  only  slightly  greater  computational  cost;  (2) 
MPI  typically  required  more  iterations  but  less  work  than  the  NILS,  due  to  the  higher 
per  iteration  work  required  for  the  NILS  solver;  and  (3)  the  IT2  transform  compared 
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favorably  with  all  other  transforms  in  terms  of  efficiency  and  robustness.  Based  upon 
these  results,  we  limited  our  additional  phases  of  this  work  to  the  KINT  and  MPI 
approaches,  although  all  transforms  were  considered. 


5.2.  Performance  Comparisons 


Using  the  general  observations  made  from  the  baseline  comparisons,  a  set  of 
simulations  was  performed  using  the  KINT  interblock  conductivity  estimate  and  the 
MPI  nonlinear  solver.  For  this  set  of  simulations,  work  and  error  were  observed  as  a 
function  of  / 3  for  all  four  transforms  and  all  eight  test  problems.  Table  4  shows  the 
allowable  transform  parameter  ranges  over  which  convergence  was  achieved,  the  optimal 
parameter  values,  and  resulting  optimal  work  and  error  values.  IT1  is  not  listed  as  a 
separate  transform  since  it  is  represented  by  the  (3  —  0  point  of  the  IT2  transform. 
Though  several  of  the  optimal  [i  values  for  IT2  are  zero  or  close  to  zero,  particularly  for 
problems  involving  nv  <  2,  significant  differences  in  performance  are  seen  for  IT1  and 
IT2  over  the  entire  range  of  problems  considered. 

The  error  vs.  work  results  are  illustrated  for  two  of  the  test  problems  —  A  and  C 
— in  Figures  3  and  4,  respectively.  The  lines  on  the  work-error  plots  represent  the  work 
and  error  for  each  transform  as  the  transform  parameter  is  varied.  Results  are  plotted 
for  parameter  values  at  equally  spaced  intervals  over  the  range  of  allowable  values.  For 
THT,  HST,  and  RFT,  data  points  are  plotted  at  equally  spaced  intervals  in  [J,  and  for 
IT2,  data  points  are  plotted  at  equally  spaced  intervals  in  log  (3. 

Based  upon  this  set  of  simulations,  the  following  observations  are  made: 


Table  4 


Figures  3  an 


1.  transforms  generally  lead  to  more  efficient  solutions,  sometimes  as  much  as  an 
order  of  magnitude  more  efficient  as  measured  in  terms  of  the  objective  function 
minimized  in  this  work; 


2.  transforms  tended  to  be  increasingly  efficient  compared  to  untransformed  solutions 


18 


as  the  sharpness  of  the  front  increased,  such  as  for  cases  in  which  a  saturated  zone 
developed  and  nv  was  relatively  large; 

3.  all  transforms  were  generally  sensitive  to  (i  for  a  given  problem. 

4.  the  optimal  value  of  / 3  varied  among  problems  for  all  transforms; 

5.  while  variations  in  efficiencies  existed  among  transforms  for  each  problem,  IT2  was 
in  general  the  most  efficient  transform;  and 

6.  IT2  was  typically  the  least  sensitive  to  changes  in  (3. 

These  data  show  that  for  sharp-front  problems,  transformation  methods  can 
significantly  reduce  computational  effort  and  increase  robustness  of  the  solution  scheme 
given  an  appropriate  value  for  (3.  However,  guidance  does  not  yet  exist  to  choose  an 
appropriate  value  of  / 3  for  any  transform  a  priori.  Further,  the  sensitivity  of  transforms 
to  media  properties  suggest  that  the  relative  efficiency  and  robustness  of  transform 
methods  will  likely  decrease  as  the  degree  of  media  heterogeneity  increases. 


5.3.  Parameter  Sensitivity  to  Discretization 


We  performed  a  discretization  study  to  investigate  the  effects  of  spatial  and 
temporal  discretization  on  the  behavior  of  all  transformation  methods  for  Problems  A, 
B,  C,  and  H  and  report  these  results  in  Table  5.  For  each  problem,  a  coarse  spatial  grid 
and  two  coarse  time-step  sizes  were  tested.  Results  from  these  numerical  experiments 
showed  that  IT2  and  RFT  were  able  to  converge  and  give  accurate  solutions  on  all  of 
the  coarse  spatial  and  temporal  grids.  For  the  untransformed  RE,  the  nonlinear  solver 
failed  to  converge  for  all  of  the  problems  except  C.  The  number  of  spatial  nodes  in  the 
coarse-grid  simulations  was  roughly  10  times  less  than  in  the  previous  simulations,  which 
resulted  in  significant  computational  savings.  Work  measure  were  scaled  linearly  as  a 
function  of  the  number  of  spatial  nodes.  Therefore,  the  IT2  and  RFT  transformations 


Table  5 
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are  able  to  provide  accurate  solutions  at  much  larger  discretization  scales,  resulting  in 
very  efficient  simulations  that  would  not  be  possible  using  the  untransformed  RE  or 
other  transformations. 

Figure  5  shows  the  solution  profiles  in  untransformed  space  for  Problem  A, 
comparing  dense  grid  and  coarse  grid  solutions  at  various  simulation  times.  The  dense 
grid  solutions,  shown  as  solid  lines,  are  achieved  using  25601  spatial  nodes  at  a  fixed 
time-step  size  of  1.56e-6  days.  The  coarse  grid  solutions,  shown  as  dashed  lines,  are 
achieved  using  41  spatial  nodes  at  a  fixed  time-step  size  of  2.0e-4.  While  the  coarse  grid 
solution  is  not  highly  accurate,  it  is  impressive  that  a  solution  was  attained  at  all  and 
this  quality  of  result  may  be  adequate  for  some  uses. 

Based  upon  these  results,  we  find  that  the  optimal  f3  is  sensitive  to  spatial  and 
temporal  discretization,  and  that  the  optimal  f3  for  a  given  discretization  may  lie  outside 
the  range  for  which  convergence  can  be  achieved  at  another  discretization,  especially 
for  the  RFT  approach.  The  general  trend  noted  was  that  the  range  of  / 3  for  which 
convergent  results  were  achieved  tended  to  become  narrower  as  the  discretization  became 
coarser.  We  also  found  that  the  objective  function  (16)  in  the  parameter  optimization 
exhibits  a  more  clearly  defined  minimum  as  the  discretization  is  coarsened,  yet  the  the 
range  of  (3  values  yielding  efficient  solutions  also  becomes  narrower. 

To  illustrate  the  sensitivity  of  optimal  / 3  to  spatial  and  temporal  discretization,  RFT 
and  IT2  [J  values  were  optimized  over  a  range  of  spatial  and  temporal  discretizations 
for  Problem  A.  These  results  are  shown  as  contour  plots  in  Figures  6  and  7.  Optimal  (3 
values  as  a  function  of  A2:  and  At  are  clearly  smoother  for  IT2  than  for  RFT.  We  also 
noted  that  for  RFT,  /3’s  closer  to  the  optimal  [i  are  required  for  good  performance.  For 
IT2,  performance  is  not  as  sensitive  to  small  deviations  from  the  optimal  (3. 


Figure  5 


Figures  6  an 
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5.4.  Parameter  Estimation 

Given  the  relative  insensitivity  of  the  IT2  parameter  to  changes  in  media 
conditions,  it  is  desirable  to  define  an  estimator  for  this  parameter  in  terms  of  known 
system  properties.  Such  an  estimator  would  be  beneficial  in  providing  an  effective 
transform  without  the  need  for  expensive  nonlinear  parameter  optimization.  Using 
finer  discretization  scales,  trends  were  detected  in  optimal  IT2  parameter  values  as  a 
function  of  media  parameters  nv  and  av.  However,  as  discretization  parameters  are 
coarsened,  the  optimal  work  and  error  values  become  more  sensitive  to  changes  in  the 
IT2  parameter.  Thus,  it  is  difficult  to  define  an  estimator  for  the  IT2  parameter  that 
will  yield  optimal  or  near-optimal  work  and  error  results  over  a  wide  range  of  media  and 
discretization  parameters  based  upon  the  set  of  results  generated  to  date. 

Preliminary  work  shows  that  it  may  be  much  easier  to  define  an  estimator  for 
the  IT2  parameter  that  will  yield  optimal  results  when  using  a  variable  time-step 
solution  method  than  for  the  fixed-time-step  cases  considered  in  this  work.  The  further 
development  of  this  notion,  or  alternative  approaches  to  estimating  a  near-optimal 
(3  as  a  function  of  media  properties,  auxiliary  conditions,  and  spatial  and  temporal 
discretization  levels  is  the  subject  of  future  work. 

6.  Conclusions 

Several  conclusions  can  be  drawn  from  our  comparison  of  transformation  approaches 
and  traditional  approaches  to  solving  RE  for  a  range  of  media  properties  and  auxiliary 
conditions. 

•  Transformation  approaches  have  the  potential  to  lead  to  more  efficient  and  robust 
solutions  of  Richards’  equation. 

•  The  potential  advantage  of  transformation  approaches  increases  as  the  difficulty  of 
the  problem  increases,  as  measured  by  the  average  number  of  nonlinear  iterations 
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per  time  step  that  are  required  for  the  untransformed  case. 

•  Capitalizing  on  the  potential  advantages  of  transform  methods  requires 
specification  of  an  appropriate  transformation  and  knowledge  of  a  reasonable  value 
for  any  free  parameters  that  exist  in  the  transform. 

•  A  set  of  one  new  (IT2)  and  four  existing  transforms  are  described  within 
a  common  framework,  and  the  optimal  transform  parameter  is  shown  to 
depend  upon  auxiliary  conditions,  media  properties,  and  spatial  and  temporal 
discretization. 

•  IT2  showed  good  efficiency  and  robustness  compared  to  all  other  transforms  for  a 
wide  range  of  media  properties  and  discretization  levels. 
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Figure  1.  Solution  profile  in  ip  for  Problem  A. 

Figure  2.  Solution  profile  in  p  for  Problem  A. 

Figure  3.  Error  vs.  work  for  Problem  A  using  MPI  solver  and  KINT  conductivity 
estimation. 

Figure  4.  Error  vs.  work  for  Problem  C  using  MPI  solver  and  KINT  conductivity 
estimation. 

Figure  5.  Comparison  of  solution  profiles  for  dense  and  coarse  grid  solutions  to  Problem 
A. 

Figure  6.  Contour  plot  of  optimal  RFT  parameter  (5  vs.  A2  and  At  for  Problem  A. 
Figure  7.  Contour  plot  of  optimal  IT2  parameter  (3  vs.  A z  and  At  for  Problem  A. 
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Table  1.  Transformation  Definitions 


Transformation 

p  Definition 

IT1 

1*00  KWW, 

— oo  <  'ip  <  oo 

THT 

4  [e«o  -  m\ + a 

ip  <  (3 

*P, 

tp>  (3 

HST 

smh[  V]* 

ip  <  /3 

f-p 

Pi  ’ 

ip>(3 

RFT 

l +W’ 

ip  <  0 

ip  >  0 

IT  2 

StooKW)  dip  +  /?  [«(</>)- 

A- 

IA 

o 

+  p(°)- 

ip  >  0 
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Table  2.  Test  Problem  A-D  Parameters 


Variable 

Problem  A 

Problem  B 

Problem  C 

Problem  D 

Or  (-) 

0.093 

0.279 

0.102 

0.000 

(-) 

0.301 

0.416 

0.368 

0.330 

av  (m_1) 

5.47 

6.12 

3.35 

1.43 

nv  (— ) 

4.264 

1.914 

2.000 

1.506 

Ks  (m/day) 

5.040 

2.298 

7.970 

0.067 

Ss  (m_1) 

1.0  x  10~6 

1.0  x  10-6 

0.00 

0.00 

Z  (ill) 

[0,10.0] 

[0,10.0] 

[0,0.3] 

[0,6.0] 

t  (days) 

[0,0.18] 

[0,0.265] 

[0,0.092] 

[0,7.16] 

ip o  (m) 

-z 

— z 

-0.10 

-0.3069 

ipi  (m) 

0.00 

0.00 

-0.10 

-0.3069 

1P2  (m) 

0.10 

0.10 

-0.75 

-0.07 

(m) 

0.0125 

0.0125 

0.0025 

0.0125 

At  (days) 

5.0  x  10~5 

1.0  x  10~4 

5.0  x  10~4 

7.16  x  10~2 
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Table  3.  Test  Problem  EH  Parameters 


Variable 

Problem  E 

Problem  F 

Problem  G 

Problem  H 

Or  (-) 

0.045 

0.057 

0.078 

0.095 

(-) 

0.430 

0.410 

0.430 

0.410 

av  (m_1) 

14.5 

12.4 

3.60 

1.90 

nv  (— ) 

2.680 

2.280 

1.560 

1.310 

Ks  (m/day) 

7.128 

3.502 

0.250 

0.062 

Ss  (m_1) 

1.0  x  10~6 

1.0  x  10~6 

0.00 

0.00 

Z  (ill) 

[0,10.0] 

[0,10.0] 

[0,5.0] 

[0,2.0] 

t  (days) 

[0,0.24] 

[0,0.45] 

[0,2.25] 

[0,1.0] 

ip o  (m) 

-z 

— z 

-z 

-z 

ipi  (m) 

0.0 

0.0 

0.0 

0.0 

1P2  (m) 

0.10 

0.10 

0.10 

0.10 

A;?  (m) 

0.0125 

0.0125 

0.0125 

0.00625 

At  (days) 

1.0  x  10~4 

1.5  x  10~4 

3.0  x  10~3 

2.0  x  10~3 

31 


Table  4.  Simulation  Results 


Problem 

Transform 

Pmin 

Pmax 

Optimal  Performance  Values 

Li  Error  Work  x 

Popt  (xlO-3)  L1  Error 

A 

UNT 

5.30 

129.85 

THT 

^Pmin 

-0.05 

-0.05 

3.93 

134.80 

HST 

'Pmin 

-0.15 

-0.15 

3.08 

97.02 

RFT 

-3.50 

0.00 

-3.50 

4.36 

140.39 

IT2 

0.00 

100.00 

0.58 

0.74 

19.24 

B 

UNT 

0.69 

17.46 

THT 

'Pmin 

-0.015 

-1.20 

0.46 

15.90 

HST 

'Pmin 

-0.04 

-1.50 

0.46 

15.89 

RFT 

-36.00 

0.00 

0.00 

0.69 

22.63 

IT2 

0.00 

20.00 

0.001 

0.50 

15.99 
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Table  4.  (continued) 


Problem 

Transform 

Pmin 

Putax 

Optimal  Performance  Values 

Li  Error  Work  x 

Popt  (xl0~3)  Li  Error 

C 

UNT 

12.14 

21.25 

THT 

Ipmin 

0.00 

0.00 

5.78 

8.84 

HST 

Putin 

-0.001 

-1.60 

6.51 

9.83 

RFT 

-10.00 

0.00 

-10.00 

5.76 

8.87 

IT2 

0.00 

100.00 

0.005 

6.27 

6.27 

D 

UNT 

0.08 

0.07 

THT 

Putin 

0.00 

-0.13 

0.08 

0.08 

HST 

Putin 

-0.001 

-0.06 

0.08 

0.07 

RFT 

-45.00 

0.00 

-45.00 

0.08 

0.08 

IT2 

0.00 

100.00 

0.04 

0.08 

0.08 
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Table  4.  (continued) 


Problem 

Transform 

fimin 

fimax 

Optimal  Performance  Values 

Li  Error  Work  x 

/3opt  (xl0~3)  Li  Error 

E 

UNT 

5.22 

110.66 

THT 

^Pmin 

-0.03 

-0.03 

0.60 

21.78 

HST 

'Pmin 

-0.05 

-0.05 

0.33 

9.87 

RFT 

-71.00 

0.00 

-71.00 

0.33 

8.51 

IT2 

0.00 

5.00 

0.05 

0.27 

7.69 

F 

UNT 

4.86 

125.87 

THT 

‘Pmin 

-0.04 

-0.04 

0.34 

13.77 

HST 

‘Pmin 

-0.02 

-0.02 

0.70 

26.46 

RFT 

-63.00 

0.00 

-63.00 

0.69 

21.25 

IT2 

0.00 

1.30 

0.01 

0.31 

11.07 
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Table  4.  (continued) 


Problem 

Transform 

fimin 

fimax 

Optimal  Performance  Values 

Li  Error  Work  x 

Popt  (xl0~3)  Li  Error 

G 

UNT 

1.61 

19.96 

THT 

^Pmin 

-0.01 

-0.90 

1.42 

23.86 

HST 

Ipmin 

-0.12 

Ipmin 

1.60 

25.76 

RFT 

-68.09 

0.00 

-47.50 

1.36 

10.74 

IT2 

0.00 

1.27 

0.002 

1.46 

14.89 

H 

UNT 

0.51 

4.04 

THT 

'Ipmin 

-0.007 

-0.01 

0.46 

5.11 

HST 

'Ipmin 

-0.22 

Ipmin 

0.51 

5.30 

RFT 

-201.00 

0.00 

-17.50 

0.45 

4.20 

IT2 

0.00 

0.73 

0.00 

0.51 

2.47 

Table  5.  Coarse  Grid  Simulation  Results 
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Media  A;? 

A  0.1 


Optimal  Performance  Values 


At 

Transform 

fimin 

fimax 

fiopt 

Li  Error 

(xlO-3) 

Work  x 

Li  Error 

le-4 

UNT 

DNC 

THT 

-2.90 

-0.50 

-0.50 

2.68 

4.50 

HST 

DNC 

RFT 

-5.49 

-2.10 

-5.49 

43.98 

68.70 

IT2 

0.00 

1.40 

0.003 

2.22 

2.95 

2e-4 

UNT 

DNC 

THT 

-0.85 

-0.31 

-0.31 

2.41 

2.73 

HST 

-1.04 

-0.92 

-0.92 

7.93 

7.88 

RFT 

-4.00 

-0.53 

-0.53 

6.08 

5.96 

IT2 

0.00 

1.80 

0.01 

2.38 

1.97 
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Table  5.  (continued) 

Optimal  Performance  Values 
Li  Error  Work  x 

Media  A;?  At  Transform  /3min  /3max  (3opt  (xlCT3)  Lx  Error 

B  0.1  2e-4  UNT  DNC 

THT  DNC 

HST  DNC 


RFT 

-71.00 

-6.80 

-50.90 

3.95 

5.17 

IT2 

0.00 

1.40 

0.00 

15.30 

30.97 

UNT 

DNC 

THT 

DNC 

HST 

DNC 

RFT 

-45.80 

-6.70 

-43.80 

3.97 

2.29 

IT2 

0.00 

1.40 

0.00 

12.68 

10.61 

Table  5.  (continued) 
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Media  A2;  A  t 

Transform 

Pmin 

Pmax 

Optimal  Performance  Values 

Li  Error  Work  x 

Popt  (xl0~3)  Ll  Error 

C  0.03  le-3 

UNT 

267.79 

9.59 

THT 

iftmin 

0.00 

-2.16 

9.89 

0.46 

HST 

^Pmin 

-0.001 

-2.51 

10.48 

0.49 

RFT 

-200.00 

0.00 

-0.28 

17.19 

0.79 

IT2 

0.00 

93.00 

0.006 

12.81 

0.45 

2e-3 

UNT 

277.55 

6.55 

THT 

'Pmin 

0.00 

-1.97 

10.71 

0.32 

HST 

‘Pmin 

-0.001 

-2.26 

12.33 

0.37 

RFT 

-200.00 

0.00 

-0.35 

18.49 

0.54 

IT2 

0.00 

93.00 

0.008 

13.61 

0.29 
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Table  5.  (continued) 


Optimal  Performance  Values 
Li  Error  Work  x 

Media  A;?  At  Transform  /3min  j3max  /3opt  (xlCT3)  Zq  Error 


H  0.05  4e-3 


8e-3 


UNT 

DNC 

THT 

DNC 

HST 

DNC 

RFT 

-33.40 

-7.03 

-17.20 

3.39 

1.35 

IT2 

0.00 

0.16 

0.00 

9.70 

2.50 

UNT 

DNC 

THT 

DNC 

HST 

DNC 

RFT 

-26.80 

-6.46 

-17.70 

3.57 

0.70 

IT2 

0.00 

0.17 

0.00 

9.48 

1.22 

•  • 


z  (m) 


Error 


•• 


Ill 


